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Abstract. We observed the pre-stellar core L1521F in dust emission at 1.2mm and in two transitions each of 
N 2 H+, N 2 D+, C ls O and C 17 in order to increase the sample of well studied centrally concentrated and chemically 
evolved starless cores, likely on the verge of star formation, and to determine the initial conditions for low-mass 
star formation in the Taurus Molecular Cloud. The dust observation allows us to infer the density structure of the 
core and together with measurements of CO isotopomers gives us the CO depletion. N 2 H + and N 2 D + lines are 
good tracers of the dust continuum and thus they give kinematic information on the core nucleus. We derived in 
this object a molecular hydrogen number density n(H 2 )~ 10 6 cm -3 and a CO depletion factor, integrated along 
the line of sight, /d = 9.5x10" /i bs(CO) ~ 15 in the central 20", similar to the pre-stellar core L1544. However, 
the A r (N 2 D + )/A'(N 2 H + ) column density ratio is ~ 0.1, a factor of about 2 lower than that found in L1544. The 
' observed relation between the deuterium fractionation and the integrated CO depletion factor across the core 

can be reproduced by chemical models if N 2 H + is slightly (factor of ~2 in fractional abundance) depleted in the 
central 3000 AU. The N 2 H+ and N 2 D+ linewidths in the core center are ~ 0.3 km s 1 , significantly larger than in 
'_h ' other more quiescent Taurus starless cores but similar to those observed in the center of L1544. The kinematical 

behaviour of L1521F is more complex than seen in L1544, and a model of contraction due to ambipolar diffusion 
is only marginally consistent with the present data. Other velocity fields, perhaps produced by accretion of the 
surrounding material onto the core and/or unresolved substructure, are present. Both chemical and kinematical 
analyses suggest that L1521F is less evolved than L1544, but, in analogy with L1544, it is approaching the 
"critical" state. 
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1. Introduction 

The pre-stellar core L1544 has recently been the subject 
of much study because while apparently in hydrostatic 
equilibrium, there are indications that it is close to the 
"critical" state at which it will become gravitationally un- 
stable and from which it will dynamically collapse (see 
discussions by Tafalla et al. ROM Ciole k & Basu I2UUU1 
Caselli et al. I^OO^al Caselli et al. I2002bjl . If correct, this 
is of fundamental importance because it defines the initial 
conditions for the formation of a protostar, which affects 
many theoretical studies of low mass star formation. 
Clearly trying to extrapolate general trends from a single 
object is difficult and a larger number of L1544-like cores 
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(preferably with the same external environment) should 
be studied. Unfortunately there are rather few other ob- 
jects with similar properties due to the short timescale of 
this phase. According to the Ciolek & Basu {212101) model, 
for example, (contraction of a disk driven by ambipolar 
diffusion) L1544-like properties fit the model structure of 
a core at times 3-8 x 10 4 years prior to the collapse, after 
an evolution of 2.6 x 10 6 years. Thus in that particular 
model, L1544 finds itself in the last few percent of its evo- 
lution prior to becoming a protostar. While this may be a 
somewhat too literal interpretation of the model results, it 
shows that "L1544-type cores" should be relatively rare. 

Further progress requires the definition of what is a 
L1544-like core. One answer is to use current estimates 
of dust emission and absorption selecting cores of dust 
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extinction upwards of 50 mag. Another approach is to 
say that cores which show signs of infalling gas (as does 
L1544, see Williams et al. H^l Tafalla et al lTMSjl are 
"L 1544- twins" . This latter indicator is complicated by 
the fact that at the high densities found in the nuclei of 
cores similar to LI 544, many molecular species and in 
particular CO and CS freeze-out onto dust grain surfaces 
(see Kramer et al. 119991 Caselli et al. 119991 Bacmann et 
al. E0I32 Bergin et al.EEEl J0rgensen et al.|20021 Tafalla 
et al. 2002}; observing such tracers implies observing 
the low density surrounding envelope. However, recent 
studies indicate that species whose abundance is linked 
to that of molecular nitrogen such as N2H + and NH3 (as 
well as their deuterated counterparts) do not condense 
out in the same fashion and hence can be used as tracers 
of the dense gas (Bergin & Langer I1997fl . The extent to 
which this is true is debatable but it is a useful hypothesis 
and substantiated by the general similarity of the spatial 
distributions seen for example in dust emission and in 
maps of N 2 H+ (Tafalla et al. 12003 IMHjl . 



Caselli et al. lP02al2OO2b|) have used N 2 H+ and N 2 D+ 
to derive the physical, chemical and kinematical properties 
of L1544. They found that L1544 has a central N 2 H+ col- 
umn density of 1.5 x 10 13 cm -2 and a column density 
ratio N(N 2 D+)/AT(N 2 H+) f 0.24. The N 2 H+ linewidths 
towards the nucleus (the dust emission peak) are roughly 
0.3 km s _1 and decrease as one goes to positions away from 
the center. The line of sight velocity measured in N 2 H + (1- 
0) and N 2 D + (2-l) shows a gradient along the minor axis 
of the elliptical structure seen in 1.3mm dust emission but 
no clear gradient along the major axis. 

In this paper, we will study another core in the Taurus 
complex, L1521F (at an assumed distance of 140 par- 
sec) using the same approach as in our study of L1544. 
Repeating the L1544 study carried out by Caselli et al. 
(2002a, 2002b) is important because it allows us to check 
to what extent LI 544 is an exceptional case. In order to 
do this we need another source which has the same gen- 
eral characteristics as L1544. The source selection was 
made using some preliminary results we obtained in a 
survey carried out at the IRAM-30m telescope. L1521F 
stood out as being the only core in Taurus, besides L1544, 
with strong N 2 D + (2-l) emission compared to N 2 H + (l-0). 
This suggests enhanced deuterium fractionation imply- 
ing an advanced evolutionary state (Caselli et al. [2002b). 
Previous observations of this object have been carried out 
by Mizuno et al. ifllMj) . Onishi et al. 1)1996[> . Codella et 
al. ifTMTjl . and Lee et al. I|1999a|l . Onishi et al. lfT^9"9"jl 
also studied L1521F (which they call MC 27), and found 
a high central density, suggesting that this is the most 
evolved starless condensation in Taurus. L1521F was also 
noted by Lee et al. p999b| as a strong infall candidate, 
in their survey of CS and N 2 H + lines in starless cores, 
although later mapping of the two tracers has shown ex- 
tended "red" asymmetry in the CS(2-1) profiles (Lee et 
al.EOOH)- 



In section 2 of this paper, we describe our observa- 
tional procedure. In section 3 we present the observational 
results deriving the physical characteristics of the source 
and analysing its chemical and kinematical properties. In 
section 4 we discuss the observational results and the sum- 
mary can be found in Section 5. 

2. Observations 

The observations were carried out between April 2002 and 
January 2003 at the IRAM-30m in three different runs. 

In April 2002, we observed the core in N 2 H + (l-0), 
N 2 H+(3-2), N 2 D+(2-l) and N 2 D+(3-2). In general, we 
used the symmetric frequency switch mode and the facility 
autocorrelator; in Tabled we summarize the main obser- 
vational parameters. The frequencies of the N 2 H + (l-0), 
N 2 D+(2-l) and N 2 D+(3-2) have been updated following 
the recent determinations of Dore et al. I|2004(l : the val- 
ues in the table refer to the F 1 F = 2 3 -> 1 2, 2 3 -> 
1 2, and 4 5-^34 hyperfine components of the N 2 H + (1- 
0), N 2 D+(2-l) and N 2 D+(3-2) transitions, respectively. 
For N 2 H + (3-2) we used the frequency of the 2 1— >1 
component as determined by Caselli et al. (2002a). In the 
case of the April 2002 N 2 D + (3-2) observations, in order to 
improve the baseline quality, we used also the "Wobbler 
switching" mode with a 240" throw. We reached an r.m.s. 
sensitivity in main-beam brightness units of about 100 mK 
in all lines except N 2 H+(3-2) (~ 400 mK). The pointing 
was checked every 2 hours by means of a 3 or 2 mm contin- 
uum scan on nearby quasars and was accurate to within 
~ 4". 

In order to refine the maps, originally taken with a 20" 
spacing, we observed in Nov. 2002 with a 10" grid (but 5" 
spacing in the inner 20 seconds). 

Table 1. Telescope settings and parameters. 



line 


frequency 


HPBW 


F throw 








GHz 




kHz 


K 


km s _1 


N 2 H+(l-0) 


93.1737725 


26 


7.5 


150 


0.063 


N 2 H+(3-2) 


279.511385 


9 


14.3 


1900 


0.021 


N 2 D+(2-l) 


154.217137 


16 


7.5 


370 


0.038 


N 2 D+(3-2) 


231.321966 


10 


14.3 


750 


0.050 


C 18 O(l-0) 


109.782160 


22 




170 


0.026 


C ls O(2-l) 


219.560319 


11 




450 


0.033 


C 17 O(l-0) 


112.358988 


21 


7.5 


230 


0.026 


C 17 0(2-l) 


224.714368 


11 


14.3 


730 


0.052 


Col. 2 line rest freq., Col. 3 Half Power Beam Width 



Col.4 Freq.Throw, Col. 5 Syst. 
Col. 6 Channel Spacing 



Temperature (K) 



Between January 2003 and March 2003, we obtained 
continuum data at 1.2mm together with observations of 
C 18 and C 17 O(l-0) and (2-1). These data were taken in 
service mode by the IRAM staff. 

The continuum data were obtained using MAMBO 
II, the 117-channels bolometer available at the 30m. We 
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mapped the core within an area of 150" x 150" scanning 
in azimuth with a 5"/sec speed and an interval between 
the subscans of 8". The atmospheric attenuation was mea- 
sured to be 0.14 based on tipping curves measured after 
the map. The data were calibrated using the sources HL 
Tau and L1551 for which we assumed fluxes of 0.9 Jy and 
1.4 Jy respectively and the final sensitivity was 5 mJy per 
10.5" beam. The calibration error inherent in this com- 
parison is likely to be at least ten percent due to both 
atmospheric fluctuations and calibration errors. 

The C 18 data were taken using the on-the-fly tech- 
nique. We simultaneously mapped the C 18 O(l-0) and 
C 18 0(2-l) using both polarizations for each line. The area 
covered was 150"xl50" and was scanned in the Right 
Ascension direction; the distance between the subscans 
was 5" as was the angular separation between two succes- 
sive dumps. We also obtained a 9 points map in C 17 O(l-0) 
and C 17 0(2-l) centered at the dust peak and spaced by 
20". 

3. Results 

3.1. Integrated intensity maps and continuum emission 

We show the observed map of the 1.2mm dust continuum 
in Fig. ^ and the map of integrated line intensity, ob- 
tained in the observed transitions of N2H + , N2D + , and 
C 18 0, in Fig. [3 The reference position for these maps is 
(04:28:39.8, 26:51:35) in J2000 coordinates. 

The bolometer map was reduced with the IRAM stan- 
dard reduction program NIC. From Fig.^ we see that the 
observed emission has a "cometary" structure in the sense 
that the low level contours are well-fitted by an ellipse 
with the maximum offset from the center. We can fit the 
general elliptical structure of the core with a 2D-gaussian 
centered at offset (-30", 20") with full width half-power 
dimensions of 274"x 170"and position angle 25°. If the 
more extended emission is not included in the fit, the 2D- 
gaussian is centered on the dust peak position at (-10", 0"), 
has half-power dimensions of 127" x 77", a position angle 
of 25° and an aspect ratio equal to 1.6. The peak intensity 
is 90 mJy/beam. 

The spectral line data were reduced using the standard 
IRAM package CLASS. A summary of line parameters at 
the dust peak is given in Table and the corresponding 
spectra are shown in Fig. |3| The spectra shown in Fig. [3] 
(as well as the values in Table [2J have been derived from 
data gaussian-smoothed to a resolution of 26" in the cases 
of N2H" 1 " and N2D+ , but unsmoothed in the case of C 18 0. 
These are also the effective resolutions of the N2H 4 " maps 
shown in Fig. [21 (the two N2D + maps are smoothed to 
16", the angular resolution at the (2-1) frequency). 

One clear result from Fig. |3|is that in L1521F there 
is evidence for an "infall signature" although less marked 
than in L1544, where it is attributed to extended infall 
onto the core (Williams et al. 1999) or to central infall in 
a depleted core nucleus (Caselli et al. 2002a). In fact, we 
see evidence for asymmetric profiles, with the blue peak 



brighter than the red peak, in N2H+ (1-0) (where we de- 
rive an optical depth of order 4 in the main component 
based on our fit to the hyperfine satellites), but the two 
peaks are not clearly separated as in the case of L1544. 

It is interesting to note that two peaks are present 
in the spectrum of N2H + (3-2) line toward the offset (- 
10,10), as shown in Fig. The hyperfine structure (hfs) 
fit to the N2H + (3-2) line, assuming the presence of two 
velocity components along the line of sight, gives Vlsr 1 
= 6.35±0.1 km s" 1 and V LSR2 = 6.55±0.04 km s _1 . The 
line widths of the two components are Avi = 0.15±0.02 
km s _1 and Av2 — 0.19±0.06 km s _1 , marginally (factor 

1.5) broader than the thermal N2H + line width at 10 
K. Although this result should be confirmed with higher 
sensitivity data, the velocities are consistent with those 
of the Nl and N2 (N2H + ) clumps observed by Shinnaga 
et al. H2004f) with interferometric observations. The hfs fit 
to the N2H + (l-0) line, assuming two velocity components 
and fixing the velocities to the values found with the (3-2) 
line, is consistent with the observed spectrum (if the two 
velocities are not fixed, the hfs fitting procedure applied 
to the N2H + (l-0) spectrum does not converge, probably 
because the two components are not well separated) . This 
suggests that the two clumps observed by Shinnaga et 
al. H2004J) in N2H + (l-0) are also present in our single dish 
data, although the N 2 H + (l-0) extended emission partially 
hides the features arising from the central region. 

Fig. |21 shows that L1521F has in common with sev- 
eral other cores (see also Tafalla et al. 120021 Caselli et 
al. I2002al2002b|l the property that while the maps of the 
nitrogen bearing molecules have a similar appearance to 
those observed in dust emission, maps in rare CO iso- 
topomers such as C 18 show no correlation with the dust. 
In fact, the distribution of C 18 O(l-0) and (2-1) is es- 
sentially flat over the area we have mapped, suggesting 
perhaps that the layer sampled in C 18 is either in the 
background or foreground relative to that seen in dust 
emission. 

We have been able to make some estimate of the op- 
tical depths of the observed C 18 lines by comparison 
with our 9-point C 17 map (20" spacing centred on (- 
5,-5)). From the integrated intensity ratio of the C 18 
and C 17 O(l-0) lines and assuming an intrinsic abundance 
ratio [C 18 0]/[C 17 0] of 3.65 (fr om [ 1 8 Q]/[ 17 Q] = 3.65, 
Penzias 1981 see Kramer et al. 1999 for a discussion of 
the validity of these techniques), we can put in general 
upper limits on the optical depth of C 18 O(l-0) of order 
unity. At the (-5,-5) and (-25,-5) offsets, there is however 
evidence that C 18 O(l-0) is optically thick and in fact we 
derive C 18 O(l-0) optical depths of 1.0 and 1.5 respec- 
tively at these two offsets (errors 50% ). In these positions 
we used the C 17 O(l-0) data and, assuming T ox = 10 K 
(consistent with the observed C 17 O(2-l)/(l-0) line ratio) 
and [ 18 0]/[ 17 0] = 3.65, we found an average C 18 column 
density of ~ 2xl0 15 cm -2 in the central 25". 

It is also interesting that we find that N2H + (l-0) 
and dust emission maps have similar sizes, whereas N2D + 
maps are somewhat smaller. Fitting a 2D gaussian to the 
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Table 2. Line parameters derived from hyperfine structure fitting at position (-10,0). The values refer to spatially 
averaged spectra (26"beam). 



line 


Vlsr 
km s _1 


AV 

km s _1 


T 


K km s" 1 


T ex 
K 


N 2 H+(l-0) 


6.472±0.001 


0.299±0.001 


17.9±0.1 


5.859±0.016 


4.94±0.03 


N 2 H+(3-2) 


6.388±0.010 a 


0.290±0.030 


1.4±1.0 


0.692±0.075 


5.1±3.7 


N 2 D+(2-l) 


6.505±0.004 


0.268±0.010 


2.2±0.4 


0.978±0.026 


4.6±0.9 


N 2 D+(3-2) 


6.507±0.006 


0.222±0.024 


1.2±0.9 


0.332±0.029 

/ onno,. 


4.1±3.2 



whereas for the other lines we used the updated values determined by 
Dore et al. (12U04I . b The integration includes all hyperfmes. 



N2H+ and N2D+ maps in the same fashion as for the dust 
emission, we find the parameters given in Table |3 for the 
angular sizes of N2H" 1 " and N2D + in L1521F. We deduce 
from these data a linear size for L1521F seen in N 2 H + (1- 
0) of roughly 14000 AU and dimensions in N 2 D+(2-l) of 
13000 x 6000 AU. It is interesting to note the smaller sizes 
derived in the higher J transitions of both species, suggest- 
ing that they are sampling somewhat higher density gas. 
The 1.2mm dust continuum on the other hand is more ex- 
tended suggesting that either an increase in excitation or 
in abundances is causing the N 2 D + and N2H + (3-2) line 
emission to increase in relative strength toward the dust 
peak. Notable also are the large aspect ratios (1.6-1.9) ob- 
served in N2D 4 " suggesting an ellipsoidal or triaxial form 
for the high density core. 

3.2. Density and mass distribution 

We expect the observed continuum emission at mm wave- 
lengths to be optically thin and can relate the 1.2mm flux 
to the H2 column density under the approximation of con- 
stant dust emissivity and temperature: 



^beam -^iv(T) /il.2mm 'Tfl 

where A^H^) is the averaged H2 column density , Si. 2mm 
is the emitted 1.2mm flux density from the core, fibeam 
is the telescope beam solid angle, Ki.2mm is the dust ab- 
sorption coefficient per gram of gas at 1.2mm, B V (T) is 
the Planck function at temperature T, and m the mean 
molecular mass. In our calculations we used a dust tem- 
perature Tdust = 10 K, close to the gas temperature found 
by Codella et al. (1997) using ammonia data (9.1 K), a 
dust opacity /ti.2mm = 0.005 cm 2 g" 1 (Andre et al. ll996J) 
and m = 2.33 amu (with uncertainties of typically 50%; 
see e.g. Bianchi et al. 2003 and Gongalves et al. I2004JI . 
This expression can be integrated over the 1.2mm map to 
derive the total gas mass. 

In this fashion, we derive a core mass from the con- 
tinuum measurements, within the large ellipse in Fig. ^ 
of 5.5±0.5 Mq corresponding to an integrated total flux 
of 3.0 ± 0.3 Jy. (The error is dominated by errors in cali- 
bration and baseline removal.) Within the 26" beam, used 
for our line measurements, the enclosed mass is O.7M . It 



is interesting to note that 5M© and O.8M are the aver- 
age total mass and 1 Jeans masss in the Taurus Molcular 
Cloud complex (Goodwin et al.]2004). 

For the purpose of comparison with model calcula- 
tions, we have used the 1.2 mm continuum data to es- 
timate the density distribution under the assumption of 
spherical symmetry. This inevitably involves a rough ap- 
proximation since the L1521F core clearly is not spheri- 
cally symmetric (the aspect ratio in the 1.2mm continuum 
emission is 1.6, see Tab. 03) and would be better approxi- 
mated with an ellipsoid. 

Nevertheless, we have followed the technique adopted 
by Tafalla et al. H2002[l and fit our data with a model of 
the form: 



^H 2 W = 







1 + 


(*)■ 







We thus circularly averaged the 1.2 mm data around the 
peak and used a % 2 routine to fit the dust profile using 
the above model integrated along the line of sight and 
convolved with a 10.5" gaussian. The result of the fitting 
procedure can be seen in Fig. El ■ The parameters that best 
fit our data are no = 10 6 cm -3 , ro = 20" (corresponding 
to 0.014 parsec or 2800 AU) and a = 2.0. 

3.3. CO freeze-out 

The comparison between the dust continuum emission 
and the C ls O integrated intensity map allows the de- 
termination of the amount of CO freeze-out onto dust 
grain surfaces, integrated along the line of sight. This 
is possible because the millimeter continuum data fur- 
nish A^(H2)i.2mm, the column density of molecular hy- 
drogen across the core, assuming optically thin condi- 
tions. The same quantity is obtained from C 18 lines 
(A^B^co)? again under the assumption of optically thin 
conditions, and adopting the relation between CO and H2 
valid in undepleted conditions ([CO]/[H 2 ] ~ 9.5xl0~ 5 = 
x(CO) can , the "canonical" CO abundance value; Frerking 
et aL lTMSjl . From the A r (H 2 )i.2mm/^V(H 2 )co column den- 
sity ratio, the integrated CO depletion factor (/d), or the 
ratio between the canonical and the observed CO abun- 
dance (a;(CO)can/a;(CO)obs), is easily derived. 

In practice, this process requires the division of the 
1.2mm dust continuum emission (jFi.2mm[mJy/beam] = 
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line Int. Aa AS Maj. Axis Min. Axis P.A. 

Kkms- 1 " 



1.2mm 




-8 


2 


127 


77 


-25 


N 2 H+(l-0) 


5.82 


-12 


1 


110 


90 


-18 


N 2 H+(3-2) 


0.68 


-7 


1 


90 


67 


+30 


N 2 D+(2-l) 


1.02 


-9 


4 


90 


46 


-9 


N 2 D+(3-2) 


0.34 


-11 


5 


63 


38 


-33 



Col. 2 Integrated intensity at peak, Col. 3 R.A. offset of peak 
Col.4 DEC offset of peak, Col. 5 Major axis size (FWHM) 
Col. 6 Minor axis size (FWHM), Col. 7 Position angle, 
measured East of North. 



£i.2mm[iriJy]/fibeam) map by the C 18 integrated inten- 
sity (Wqisq) map. The map-division has been carried out 
using the IRAM image manipulation software GRAPHIC, 
after degrading the continuum map to the angular reso- 
lution of the C 18 O(l-0) observations (22") so that the 
integrated depletion factor can be expressed by: 



/ D = ir(CO) can x 



JV(H 2 )i 



8.5 x 10" 



A^(C 18 0) [ 16 0] 
2^1. 2mm (mJy/22"beam) 



(1) 



W^is^Kkms- 1 ) 

To derive eq. Q we have assumed a dust temperature 
Tdust = 10 K, Ki.2mm = 0.005 cm 2 g -1 , and an abundance 
ratio [ 16 0]/[ 18 0] = 560 (Wilson & Rood Q3gg| . 

There are several caveats to the above procedure. One 
is that the "canonical abundance" appears to vary from 
cloud to cloud (Lacy et al. 1991 Alves et al. 1999. Kramer 
et al. I1999[> and is roughly one third of the diffuse cloud 
carbon gas phase abundance (Sofia et al. I1997|) . Given 
that CO is expected to represent essentially all the gas 
phase carbon in molecular clouds, this suggests depletion 
of carbon in some form (not necessarily as solid CO) even 
on the outskirts of cores (we note that the direct study 
of CO solid state features in Taurus (Chiar et al. I1995f) 
shows that solid CO towards 4 field stars in Taurus is less 
than 40% of the canonical gas phase CO abundance and 
is observed for extinctions Ay greater than 6 magnitudes). 
Thus we conclude that in particular cases such as L1521F, 
it is quite possible that we are using a value of x(CO) can 
which is a factor of order 2 too large or small thus influ- 
encing the values of /d which we derive but not the trends 
over our map. 

Another problem is that the values of /d which we 
derive are integrated along the line of sight in a situa- 
tion where the observed C 18 emission forms in an outer 
shell whereas the dust emission mainly emanates from the 
dense core nucleus. As a consequence, we observe C 18 
mainly from foreground and background gas which (see 
Fig. EJ) has an essentially flat distribution over the region 
mapped by us. One concludes that our values of /d are 
strict lower limits to the CO depletion in the core nucleus 
from which dust emission is observed. It is also the case 
that in this situation, the form of our map of fjj is es- 
sentially that of the dust emission (as we indeed find, see 



eq. EJ- The local distribution of the CO depletion factor 
(which we call JFd, see below) may differ greatly and be 
much more highly peaked than in our map. Nevertheless, 
the /d values derived by us are direct observables and 
we have therefore used them for the purpose of correlat- 
ing with parameters such as the observed deuterium frac- 
tionation. We have also used them for model comparisons 
(sections 14. 1 . II and 14.1.21) bearing in mind the above dis- 
cussion. 

In Fig. the /d map is shown (thick contour) over- 
lapped with the smoothed 1.2mm map (grey scale) and 
the N 2 D+(2-l) map (dashed contours). The CO deple- 
tion factor increases between 6 at the lowest contour of 
the 1.2mm map (.Fi.2mm = 95 mJy/22" beam) to 18 at 
the peak position (offset [0,9]), which is 11" off the mil- 
limeter dust emission peak (offset [-8, +2]), where ^"l^mm 
= 320 mJy/22" beam). From Fig. we immediately see 
that /d correlates with the N2D 4 " emission (and deuterium 
fractionation, see Sect. 13. 4|) and the 1.2 mm dust flux. The 
good correlation between /d and Ji.2mm is also evident in 
Fig. where /d is plotted as a function of A^H^i. 2mm 
(= 4.25xl0 20 J 7 i.2mm mJy/(22" beam), with the assump- 
tions on Tdust and Ki.2mm as described above). Thus, in 
L1521F, and with the caveats discussed above, /d is lin- 
early dependent on the H 2 column density (-/V(H 2 )); the 
best fit to the data in Fig. (dashed curve) gives: 



/d 



1.5 



N(H 2 



10 22 cm 2 



(2) 



In retrospect, this is not surprising as N(C ls O) does not 
vary greatly over the map. We note that /d is always ^ 
2 in the region traced by the present observations, sug- 
gesting that even at the core edges (where A^H^) ~ 10 22 
cm~ 2 ) CO molecules are partially (~ 30%) frozen onto 
grain mantles. This result is consistent with the average 
value of CO depletion (~ 25%) gauged from observations 
of solid CO features in the direction of background stars 
in the Taurus complex (Chiar et al. 1995). The thick curve 
in Fig. dis the result of a simple chemical model of a cen- 
trally concentrated cloud, where CO depletion is taken 
into account, and where the CO abundance has been in- 
tegrated along the line of sight to obtain the observed /d 
(/ D = / F D (r) dl/ J all, with F D (r) being the CO depletion 
factor within the core, thus function of the cloud radius r; 
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see Section |4~T)) . We anticipate here that the /d value at 
the cloud center is only a very small fraction (~ 1%) of 
-Fb in the central few thousand AU of the core. 

3.4. N2H + - N2D + column densities, deuterium 
fractionation, and volume densities 

The N2H + and N2D + column densities have been de- 
termined using the "constant-T ex " (CTEX) approxima- 
tion, which reduces to simple analytic expressions (see 
Appendix A in Caselli et al. [2002b), and the Large 
Velocity Gradient (LVG) approximation. Both approaches 
give reasonable column density estimates as long as opti- 
cal depths are small. When, as for example for N2H + (1- 

0) , the optical depths in the main hyperfme satellites are 
large, one is best (independent of method) to use the 
weakest of the satellites or alternatively the optical depth 
inferred from the intensity ratio of the weakest satellite 
to the strong main components. The errors in any case 
stem from the difficulties in estimating the optical depth 
of thick lines compounded with possible non-LTE effects 
for the hyperfme satellites (Caselli et al. 1995 ). Errors due 
to the estimate of the partition function (i.e. the fraction 
of the species in unobserved levels) appear to be less. We in 
general report column densities for N2H + using the CTEX 
approach based on the integrated intensity of the weakest 
satellite of N2H + (l-0) and for N2D + assuming optically 
thin conditions. From comparison between the different 
approaches employed by us, we estimate the column den- 
sity errors to be 30 percent. 

We have also used LVG estimates to infer the density 
towards the peak and edges of our N2H + map. Here, we 
assume a temperature of 10 K and have used rates from 
Green {1975) for collisions of He with N2H 4 ". Based on 
the values in Tab. [21 (data smoothed to 26" spatial resolu- 
tion) for the integrated intensities of N2H + , we find n(H2) 
= 4 x 10 5 cm~ 3 towards the dust peak of L1521F and 
2.5xl0 5 cm- 3 30" North (offset [-10, 30]). For N 2 D+(2- 

1) and (3-2), the corresponding numbers are 6xl0 5 cm -3 
and 3.5xl0 5 cm -3 , at the peak and (-10, 30) offset posi- 
tions, respectively. These values are similar to the corre- 
sponding values for L 1544 consistent with the idea that 
they have similar density distributions. The density esti- 
mates are somewhat smaller than estimates based on the 
dust emission, probably due to the different (factor of 2.4 
lower) spatial resolution 1 , and consistent perhaps with the 
idea that N2H 4 " is abundant in a shell outside but close 
to the dust peak. However, given that the LVG method 
assumes homogeneous conditions, the LVG-derived densi- 
ties are averages along the line of sight, thus lower values 
are expected when compared to the continuum data anal- 
ysis, which takes into account the core density structure. 

The deuterium fractionation is directly estimated from 
the iY(N 2 D + )/V(N 2 H + ) column density ratio (= i? dcut ), 
and the i?dcut map in L1521F, assuming CTEX condi- 

1 The central density obtained with the continuum data re- 
duces to ~ 6xl0 5 cm -3 when averaged within 26". 



tions, is shown in Fig. OH -Rdeut ranges between ~ 0.02 at 
the core edge to 0.1 in a region about 20" in size and cen- 
tered on the dust peak position. The peak value of i?dcut 
is about a factor of 2 smaller than that found in L1544 
(Caselli et al. I2002b|) . We note that the column density 
of N2H + is similar in the two cores, and that the factor 
of 2 of difference in deuterium fractionation is due to the 
(factor of 2) larger N2D + column density in L1544. This 
suggests that, although the two cores are similar in struc- 
ture, L1521F is probably slightly less evolved than L1544 
(see Sect.UJl. 

3.5. Line width vs. impact parameter 

Quiescent starless cores mapped in NH3(1,1) and 
N2H + (l-0) lines typically show a "velocity coherent" 
zone of nearly constant line width (Av) within the half- 
maximum contour map, followed by a Av rise at larger dis- 
tances from core center (e.g. Goodman et al. H998ll . There 
are however exceptions to this general trend, as pointed 
out by Caselli et al. l|2002c[) . In particular, L1544 shows a 
significant increase of N2H + and N2D + (but not H 13 CO + 
and DCO + ) line widths toward the center (factor of 1.5 
within ~ 50"; Caselli et al. 12002a} . This increase has been 
interpreted as evidence of infall in the central few thou- 
sands AU, where CO and related species (such as H 13 CO + 
and DCO + ) are heavily depleted. Indeed, the line-width 
increase is consistent with models of ambipolar diffusion 
(see Sect. Ell . 

The same trend has been observed in L1521F using 
N2H + and N2D 4 " lines, and in Fig. |5|we show the results 
obtained in N2H + (l-0) lines (the linewidth corrected for 
optical depth, or the intrinsic linewidth, is plotted). The 
figure also shows two theoretical predictions which will be 
discussed in Sect. 14.21 The decrease of the intrinsic line 
width with impact parameter in L1521F, although not as 
steep as in L1544, is clear in Fig.El where the average value 
of Av within bins of 15" (see big dots) ranges between 0.3 
km s _1 at the dust peak to 0.25 km s _1 at a projected 
distance of 80". 

One might interpret this line width gradient as be- 
ing due to increased optical depth towards the dust peak. 
This however seems unlikely as illustrated in Fig. 1101 
where the profiles of the isolated hyperfme component of 
the N2H + (l-0) line along the 45° and -45° axes, passing 
through the dust peak position, are shown. If the line is 
self-absorbed toward the core center, the hyperfme com- 
ponents with the largest statistical weight will be broad- 
ened compared to the weakest ones, affecting the hfs fit. 
However, we performed gaussian fits to all the hyperfines 
components, finding the same Av values within the un- 
certainties. In fact, a similar trend is also observed if one 
plots the linewidth of the Fi F = 1 ^ 1 1 (or weakest) 
component of the N2H + (l-0) line versus the impact pa- 
rameter. This line, being thin and symmetric across the 
core, is not affected by self-absorption. 
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We believe that this is a characteristic of unstable or 
supercritical (e.g. Mouschovias & Spitzer 1976J cores on 
the verge of star formation, more briefly called pre-stellar 
cores, a term which we use to characterize the subset of 
starless cores undergoing central infall (see also Ward- 
Thompson et al. l2UU3|) . 

3.6. Velocity field 

Assuming that L1521F is in solid body rotation, we de- 
termined the magnitude Q and the direction 8 of the cor- 
responding velocity gradient following the procedure de- 
scribed in Goodman et al. ( 1993). The magnitude Q ranges 
between 0.4 and 1 km s _1 pc _1 depending on the tracer 
used, and large variations are also obtained for the gra- 
dient direction O (see Tab.0J columns 2 and 3). There is 
no tendency for Q to increase for higher density tracers, 
as observed in L1544 (Caselli et al. I2002af) . 



Table 4. Results of velocity gradient fits. 



line 


Q 


e a 


<Gi> b 




(km/s/pc) 


(deg) 


(km/s/pc) 


N 2 H + (l-0) 


0.366±0.005 


-133.6±0.7 


1.9±0.1 


N 2 H+(3-2) 


l.lzbO.l 


-168±9 


2.7±0.3 


N 2 D+(2-l) 


0.95±0.06 


+145±5 


3.5±0.2 


N 2 D+(3-2) 


0.4±0.2 


-128±27 


4.7±0.4 



NOTE: a Direction of increasing velocity, measured East 
of North. b Mean values of the magnitude of local velocity 
gradients and corresponding standard error. 



To investigate in more detail the internal motions of 
L1521F, we determined "local" velocity gradients (see 
Caselli et al. 2002iJ), where Gi and 9i have been calcu- 
lated in adjacent 3x3-point grids of the maps. The re- 
sults are shown in Figs, [ill and El for N 2 H+ and N 2 D+ 
maps, respectively. The arrows across the map indicate 
the magnitude and the direction of local gradients, and 
they are centered on the 9-point grid used to estimate 
the corresponding Q\ and Oi values. From the figures it is 
clear that L1521F is not undergoing solid body rotation. 
The velocity structure is quite complex, showing portions 
of the core where the gradient direction changes rapidly. 
For example, if we concentrate on the N2H + (l-0) map 
(see Fig. ^2 [top panel]), which has the highest sensitiv- 
ity, three converging velocity gradient patterns are clearly 
visible: (i) toward South-West in the Northern half, (ii) 
toward North-East in the SE quadrant, and (iii) toward 
North-West in the SW quadrant of the core. A similar 
structure is also present in the other maps (see Figs. 1111 
[bottom panel] and H2*| . 

Interestingly, the mean of the local-gradient mag- 
nitudes (< Q\ >) increases going from N2H + (l-0) to 
N2D + (3-2) (see column 4 of Tab. suggesting that in- 
ternal motions, although complex (0i widely varies across 
the cloud), tend to become more prominent toward the 
core center. Moreover, N2D" 1 " gradients appear larger than 



those derived from N2H" 1 ". If the observed velocity struc- 
ture is at least partially due to inward motions in the core 
center, the larger local-gradient magnitudes detected in 
N2D + can be explained by N2D + being a better tracer 
of the core central regions, as found in L1544 (Caselli et 
al. I2002al b: see also Sect. 14. l(l . However, the magnitude of 
the local gradients is generally within a few km s _1 pc -1 
(see Tab.0}, thus they are produced by small velocity vari- 
ations (Sv ~ 0.02-0.05 km s _1 ) within ~ 0.01 pc, the size 
of the grid where local gradients have been estimated (see 
also Sect. 14.21 for discussion). 

We note that the velocity gradient found with the 
present N2H + (l-0) observations is significantly smaller 
than that deduced by Shinnaga et al. ( 2004) using inter- 
fcrometric observations, probably because they resolve out 
the more extended emission, with smaller or almost oppo- 
site (see Shinnaga et al. 2004) velocity gradients, com- 
pared to the central regions. 

4. Discussion 

The previous sections described the results found in our 
analysis of L1521F. The density profile, the CO depletion 
factor, the deuterium fractionation and the velocity field 
in the core have been presented. We found several sim- 
ilarities between L1521F and L1544, including the den- 
sity structure (with n(H2) ~ 10 6 cm -3 at the center), the 
amount of integrated CO depletion (/d — 15) toward the 
dust peak, and the decrease of line width with increasing 
distance from the cloud center. Differences are however 
present in the amount of deuterium fractionation (factor 
of ~ 2 less than in L1544) and in the velocity structure. 
In this section we will discuss these findings in view of 
our knowledge of the chemical and physical processes in 
dense cloud cores. The discussion is thus split in two sub- 
sections, one focussing on the chemistry and the other on 
the kinematics of LI 52 IF. 

4.1. Chemistry 

In cold clouds, the exothermicity of the proton-deuteron 
exchange reaction 

H+ + HD -> H 2 D+ + H 2 + AE, (3) 

where AE/k = 230 K (e.g. Millar et al. IT3g9")> . is respon- 
sible for the enhancement of the H2D + /H^ abundance 
ratio well above the local elemental abundance of deu- 
terium (1.5xl0 -5 ; Oliveira et al. l2003[) . As a consequence, 
given that H2D + transfers its deuterium to neutral species 
such as CO and N2, producing DCO + and N2D" 1 ", the 
N(DCO + )/N(HCO + ) and iY(N2D+)/Ar(N 2 H+) column 
density ratios reach the values observed towards dense 
cloud cores (~ 0.02-0.2; e.g. Bu tner et al. lTMBI Williams 
et al.dSl Caselli et al. I2002bl Sect.EOjl. 

However, the freeze-out of neutral species, such as 
CO, O, and N 2 , boosts the deuterium fractionation (e.g. 
Dalgarno & Lepp[1984). In fact, a decrease in the frac- 
tional abundance of gas phase CO leads to a decrease of 
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the H3 and H2D+ destruction rates and to an increase 
(caused by the higher abundance) of the H2D+ forma- 
tion rate (e.g. Roberts & Millar 2000al see reaction[3J) • An 
empirical relation between CO depletion and deuterium 
fractionation in prestellar cores has been determined by 
Bacmann et al. (2003) using doubly deuterated formalde- 
hyde, whose abundance is also predicted to largely in- 
crease with CO freeze-out (Roberts & Millar 2000bJ see 
also Roberts et aL EfflEty . 



In the previous sections, we found that in L1521F, CO 
is depleted (with percentages ranging from 30% at the 
core edge to 93% at the center, deduced from the inte- 
grated CO depletion factor) and, similarly to L1544, the 
deuterium fractionation is large. The present estimates of 
i? deut ee iY(N 2 D + )/7V(N 2 H+) as a function of distance 
from the dust peak allow us to study the relation between 
deuterium fractionation and CO freeze-out across a dense 
core for the first time and test chemical models. In Fig. 1131 
-Rdcut is plotted as a function of /d for all the positions 
present in Fig. [S] We note a clear tendency for the deu- 
terium fractionation to increase with integrated CO de- 
pletion factor. 



In the following, we will analyse the observational re- 
sults with simple chemical models, to better understand 
the observed /d - ^!^) and -Rdeut - /b trends shown in 
Figs. 171 and 1131 and the differences between L1521F and 
L1544. It will be interesting to see whether in L1521F 
there is evidence of the so-called "molecular hole" , or the 
region where all species heavier than helium are heavily 
(*S 98%) depleted from the gas phase, as deduced in the 
case of L1544. In fact, the recent detection of H2D + to- 
ward the L1544 dust peak is consistent with the presence 
of a molecular hole in the central ~ 2800 AU (Caselli 
et al.EUUS Walmsley et al. EtM}) . Then, the N 2 H+ and 
N 2 D + emission maps should show a central hole or emis- 
sion plateau of similar size, but this has not been observed 
perhaps because of the poor spatial coverage of the cen- 
tral region and the limited spatial resolution (~ 2000 AU; 
Caselli et al. 2002a). Such a "molecular hole" was pre- 
dicted by Caselli et al. lfl"§9"9"|l and Caselli et al. I|2002afl in 
an attempt to interpret the origin of double peaked opti- 
cally thin lines for the case of L1544, including the thinnest 
N2H + (l-0) hyperfine component. In the case of L1521F, 
molecular lines, although asymmetric, do not clearly show 
two separated peaks. Thus, if the line asymmetry is due to 
the presence of a molecular hole in a contracting core, the 
hole should have a smaller size than in the case of L1544. 
Evidence for N2H + depletion toward core centers has also 
been claimed in B68 (Bergin et al EUU2) and L1512 (Lee 
et al. 12003(1 . two starless cores with central densities ~ 10 5 
cm~ 3 and close to hydrostatic equilibrium, thus in a dif- 
ferent chemical and dynamical phase compared to LI 544 
and L1521F. 



4.1.1. A simple analytical model 

The four thin curves in Fig. ^] refer to the outputs of 
a simple steady state analytical chemical model includ- 
ing H+, H 2 D+, H 2 , HD, CO, N 2 , HCO+, DCO+, N 2 H+, 
N 2 D + , electrons and negatively charged grains. The re- 
combination on grain surfaces uses the rates from Draine 
& Sutin ( 1987), assuming that the grains are bare and that 
their abundance by number is given by the MRN (Mathis 
et al. 11977(1 distribution with upper cutoff radius of 0.25 
fim and lower cutoff radius a min = 50 A(standard case, 
solid curves), and that all grains are negatively charged (a 
more realistic value for the fraction of negatively charged 
grains may be ~ 0.5; see Flower & Pineau des Forets|5003 ) . 
We also considered a larger a m i n (= 500 A, dashed curves) 
to roughly take into account the process of grain coagula- 
tion in dense cores (e.g. Ossenkopf & Hcnning 1994). As 
expected, larger grains cause a (slightly) higher deuterium 
fractionation, given that the number density of dust grains 
decreases and does also the recombination rate on grains 
of H 2 D + molecular ions (further increasing the grain size 
does not significantly change the result, given that dis- 
sociative recombination becomes more important). This 
simple model also assumes that the electron fractional 
abundance x(e) varies with density as 1.3x 10~ 5 n(H2)~ ' 5 
(McKeeHHEHl 2 , and that n(H 2 ) = 9.9 x 10 3 f^ A (in the 
density range between 3xl0 4 cm" 3 and 10 6 cm -3 ), as 
found from a linear least square fit to the observed /d 
data in Fig. and the n(H 2 ) data derived in Sect. 13.21 
from the 1.2 mm continuum observations. This allows us 
to find a simple relation for the deuterium fractionation 
as a function of /d in "L1544-like" cores, if one neglects 
multiply deuterated forms of Hjj" (see below) : 

= [N 2 D+] ^ l/3[H 2 D+]/[H+] 
^ dcut " [N 2 H+] " l + 2/3[H 2 D+]/[H+]' [ } 

and 

[H 2 D+] = fc HD [HD] 

[H+] k CO [CO] + % 2 [N 2 ] + k e [e] + k gr [gr] ' 1 j 

where fcjjD is the rate coefficient of reaction [21 (see below), 
[HD]/[H 2 ] ee 2x[D]/[H] = 3xl0~ 5 (Oliveira et al . ET)T)3Tl . 
fc co = 3.3 xl0~ 9 (T/10K)- - 5 cm 3 s" 1 is the rate coeffi- 
cient of the reaction H 2 D + + CO — > products , fcj^ 2 = 
1.7xl0~ 9 cm 3 s" 1 is the rate coefficient for the reaction 
H 2 D+ + N 2 ^ products, [N 2 ]/[H 2 ] = 2xl0~ 5 (Lee et al. 
ITMty . k e = 5.5xl0" 7 (T/IOK)" - 65 cm 3 s" 1 is the dis- 
sociative recombination rate of H 2 D + (see Caselli et al. 
1998|, and k gr is the rate of recombination of H 2 D + on 
negatively charged dust grains, which, following Draine & 
Sutin 1(1987(1 . is given by: 

McmV 1 ) = 1.6 x 10- 7 am R in ( —) 
9 V ; 10- 8 cm V10K7 

2 The use of the relation found in L1544 by Caselli et al. 
(2002b) shifts the theoretical curves up by a factor of ~ 1.3 in 
deuterium fractionation. 
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X fl + 3 .610- 4 -^--^f^V (6) 
\ 10 K 10~ 8 cm / 

Note that k gr x [gr] is equivalent to a gr in Draine & Sutin 
( 1987), where [gr] = n gr /n(H 2 ) is the fractional abundance 
of dust grains and n gr = J n gr (a)da. Introducing numerical 
values into eq.0 we obtain an expression for [BsD+J/fHg ] 
which only depends on /b and a m i n : 

[H 2 D+] = (fcHp/1.5 x lQ-Ws- 1 ) 

[H+] 7// D + 2.8 + 1.6//°- 7 + 37.8 /(a min ) ' [ ' 

where 

/(flmin) = (l.6x 10-Ws" 1 ) (lO "cm) ' (8) 

The last term in the denominator of eqn. {7\ can be ne- 
glected as soon as a m ; n = 0.05 /im or larger, so that, for 
molecular ions, dissociative recombination becomes more 
important than recombination on grains. Hence, it can be 
neglected in dense clouds, where small grains are likely to 
be depleted on the surface of large grains or coagulated in 
larger fluffy structures (Ossenkopf 1993). Depletion itself 
causes a m i n » 0.01 fim (Walsmley et al. l2004fl . 

The value of fcnD typically used in chemical codes 
is 1.5xl0~ 9 cm. 3 s , which we call the "standard rate" 
(see e.g. Roberts et al. I2003|l . However, Gerlich et al. 
( 2002J) have recently measured a lower value for this rate 
(3.5X10" 10 crmV 1 , the "GHR rate") and in Fig. [H re- 
sults obtained with the "standard" and "GHR" rate are 
shown. The observed data points lie between the two 
curves, and the best fit (dotted curve in bottom figure) 
is obtained with fcjjD — 7.5xl0 -10 cm 3 s _1 , which may 
suggest that a rate slightly (factor or ~ 2) larger than the 
one recently measured is probably needed to explain ob- 
servations. However, one should note that in the case of 
L1544 (iZdeut ~ 0-24 at /d ~ 10), even the "standard rate" 
fails to reproduce the large deuterium fractionation ob- 
served toward the dust peak. The only way to reach i?deut 
~ 0.24 with this analytical model (dash-dotted curve in 
Fig. I13JI is to allow a drop in the abundance of N 2 in the 
central regions where /d > 10 and include in the chemi- 
cal scheme D 2 H + ( 3 ), which becomes abundant in heavily 
(CO, N 2 , and O) depleted regions (Roberts et al. 120031 
Walmsley et al. 2004). Therefore, the difference in deu- 
terium fractionation between L1521F and L1544 is likely 
to be due to a different evolutionary stage, with L1521F 
being less evolved than L1544 (see also Aikawa et al. 2003). 
The inclusion of the reaction H 2 D+ + HD — ► D 2 H+ + H 2 
(Gerlich et al. I2002|) in this model leads to an increase of 
the deuterium fractionation by a factor of 2-3. However, 
in this simple chemical scheme we did not include the so- 
called "back" reactions due to ortho-H 2 (see Gerlich et al. 
2002 ), which have the effect of reducing the deuteration 

3 When D2H + is included in the model, equation 0] becomes: 

„ „ i/3[H 2 D+]/[H+]+2/3([H 2 D + ]/[H+]) 2 
dcut " i+2/3[H 2 D+]/[H+] + i/3([H 2 D+]/[H+] )2 - 



degree (see Walmsley et al. I2004fl ; this analysis is beyond 
the scope of the present paper. 

4.1.2. A simple chemical model in a centrally 
concentrated core 

The analytic calculation outlined in the previous section 
has at least two major defects. One is the assumption of no 
abundance variations within the dense core, which should 
be computed first in order to estimate molecular abun- 
dances as a function of radius and then determine the col- 
umn densities via integration along the line of sight. The 
second is our neglect of reactions with species such as O 
which also act to limit deuterium fractionation. Here, we 
present a small toy model aimed at overcoming these prob- 
lems, and already described in Casclli et al. (2002b). This 
model follows the (steady-state) chemistry of a spherically 
symmetric cloud with a density profile deduced from the 
1.2mm dust continuum emission (see Sect. l3~2f . and con- 
stant temperature T = 10 K. Neutral species in the model 
are H 2 , CO, N 2 , and atomic oxygen. We follow depletion of 
these species onto dust grains and their desorption due to 
the cosmic ray impulsive heating of the dust, following the 
procedure by Hasegawa & Herbst (1993). The abundance 
of molecular ions such as HCO + , N 2 H + and correspond- 
ing deuterated isotopomers are given by the steady state 
chemical equations using the istantaneous abundances of 
neutral species. Analogously, the electron fraction x(e) can 
be computed in terms of global estimates for the molecular 
and metallic ions and using the same istantaneous abun- 
dances of CO, N 2 , and O in the gas phase. The approach 
we have adopted is a simplified version of the reaction 
scheme of Umebayashi & Nakano (1990), which includes 
molecular ion recombination on grain surfaces using rates 
from Draine & Sutin ifTMTjl (see Caselli et al. I2002bl for 
details). 

This model furnishes the abundances of gaseous 
species as a function of radius, and the corresponding col- 
umn densities as a function of impact parameter are calcu- 
lated taking into account the appropriate beam convolu- 
tion to simulate the observations. We used the new value 
of the dissociative recombination rate for Hjj" (4x 10 -7 cm 3 
s _1 at 10 K) determined by McCall et al. I|2003|) . assumed 
fcnD = 3.5xl0~ 10 cm 3 s~ x (see previous section) and var- 
ied other parameters (essentially the binding energies on 
grain surfaces of N 2 and O). The model has been forced 
to reproduce within 10% the observed C 17 0, N 2 H + , and 
N 2 D + column densities toward the dust peak position 
(5.5xl0 14 , 1.6xl0 13 , and 1.5xl0 12 cm -2 , respectively), 
and to reproduce within a factor of 2 the observed column 
density profiles in the above molecules. The best fit bind- 
ing energies for N 2 and O (800 K and 750 K, respectively) 
are quite close to the values deduced from theoretical cal- 
culations and laboratory measurements (800 K and 750 
K, for N 2 and O, respectively; see discussion in Hasegawa 
et al. ESS and Bergin & Langer H537|l . 
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As illustrated in Fig. O the fractional abundance of 
N2H + decreases from about 2xl0 -10 at the edge of the 
cloud to about lxl0~ 10 at the center. The increase of the 
H 2 D + /H^~ abundance ratio toward the center boosts the 
formation of N 2 D + , which presents an abundance increase 
of about one order of magnitude from the edge (where 
ir(N 2 D+) = 3xl(T 12 ) to the central 4000 AU 0(N 2 D+) = 
2x 10 -11 ) of the core. Therefore, in the case of L1521F, we 
do not need to invoke any "molecular hole" as in the case 
of L1544 (see Caselli et al. EES), although the present 
data do not allow us to distinguish between the models 
with different amounts of N 2 freeze out (and consequent 
N 2 H + depletion) within few thousands AU. On the other 
hand, recent observations of ori/io-H 2 D + toward this ob- 
ject (Caselli, van der Tak, Ceccarelli et al., in preparation) 
strongly suggest that the molecular hole in LI 52 IF must 
be smaller than in L1544, given that the line is about 
two times less intense than in L1544. Indeed, we predict 
here an H 2 D + abundance of 3xl0~ 10 at radii less than 
3000 AU, a factor of about 3 lower than in L1544 (see 
Caselli et al. 2Q03J - This is another indication that L1544 
is likely to be more evolved, and more centrally concen- 
trated, than L1521F. This seems to contradict the ob- 
servational evidence that the central densities in the two 
cores are quite similar. However, one should keep in mind 
that the central density values (and the density structure) 
in the two cores have been obtained assuming isothermal 
conditions. Allowing the temperature to decrease toward 
the center, as predicted by models of dense clouds heated 
by the external radiation field (e.g. Galli et al. I2()()2|l one 
finds that lower values of the central temperature implies 
larger central densities (e.g. Evans et al. 120011 Zucconi 
et al. 1200111 . One possible solution to the L1521F/L1544 
puzzle is that L1544 is colder and more centrally concen- 
trated than LI 52 IF and hence is more depleted in the 
nucleus. This is consistent with LI 544 being more evolved 
and closer to the "critical" state than L1521F. 

We note that the chemical composition shown in 
Fig. ^] reproduces the observed Rdeut-fo (see the thick 
curve in Fig. I13|l and /d _ -^(H 2 ) (see the thick curve in 
Fig. E| relations, without any need to change the value of 
&hd (see previous section). This demonstrates the impor- 
tance of taking into account core structure and differential 
molecular freeze-out in chemical models. Moreover, note 
that the depletion factor within the cloud, Fd, is signifi- 
cantly larger (more than two orders of magnitude) than 
the integrated CO depletion factor /d, which is "diluted" 
along the line of sight (compare Fd in Fig. ^] with /d in 
Fig-0- Finally, the central value of the electron fraction 
is ~5xl0~ 9 , implying an ambipolar diffusion time scale 
(see e.g. Shu et al. HQEHl of - 2.5xl0 13 ir(e) = 1.2xl0 5 
yr, only slightly (factor of ^3) larger than the free-fall 
time scale, once again suggesting that the core is close 
to dynamical collapse (although not as close as L1544). 
As in L1544, the major ion in the core is HsO + , which 
is due to the presence of a significant fraction of gaseous 
atomic oxygen in the model (see also Aikawa et al. [2001 
for similar results). We should however point out that the 



more recent models of Aikawa et al. H2003fl . where surface 
chemistry is included, predict a much lower abundance of 
O in the gas phase, given that in this model, an O-atom 
sticking to a grain is converted to water, which remains 
on the surface (assuming desorption due to the cosmic-ray 
impulsive heating of dust grains; see Hasegawa & Herbst 
11993(1 . Low ionization potential elements, essentially S + , 
Mg+, Fe+, Si+, Na + , generically called "metals" (M+ in 
Fig. Ibfj) . are assumed to freeze out onto dust grains at 
the same rate as CO molecules. For this reason, they are 
negligible carriers of positive charges within the core, in 
our model. 

4.2. Kinematics 

In order to facilitate the comparison between L1521F and 
L1544, we analyzed the kinematical properties of L1521F 
using the same models as in Caselli et al. (2002a). In par- 
ticular, starting from the velocity field predicted by the 
ambipolar diffusion models of Ciolek & Basu ( 2000 here- 
after CB) , we have derived the linewidth gradient and the 
line profile in a disk-like geometry and compare the results 
with our observations. 

From the analysis of N 2 H + (l-0) data, we found that 
the line width, similarly to L1544, decreases with distance 
from the dust peak (see Fig. EJ). As seen in Caselli et al. 
(2002a), this observational evidence is consistent with the 
predictions of the CB model. Here, we applied the kine- 
matic analysis of L1544 (Caselli et al. I2002af) to the case 
of LI 52 IF, assuming that the cloud has a disk-like shape 
and is centrally concentrated, with the center coincident 
with the 1.2mm dust continuum map peak. From the core 
axial ratio, and following equation (1) of CB, the angle 
between the vertical axis of the model and the plane of 
the sky is found to be 18°. The disk is contracting via 
ambipolar diffusion and the resultant velocity field is used 
as input in a model which computes synthetic profiles of 
optically thin lines for all lines of sight across the model 
disk (for details, see Caselli et al. I2002af) . As in the case 
of L1544, we also assumed that the density profile and 
the radial velocity field is given by the CB model at time 
t = 2.66 Myr (= t-j in CB), which best reproduces the 
continuum observations of both cores. Line broadening is 
both due to thermal motions (uth = 0.05 km s _1 at 10 K) 
and microturbulence described by a turbulent velocity at u 
independent of position. 

Fig. 1151 shows the comparison between observed (hys- 
togram) and synthetic (curves) profiles of the N 2 H + (l-0) 
isolated hyperfine component (F±F = 01 — ^12) along the 
minor and major axes of the L1521F core. The (optically 
thin) weakest component (F±F — 10— >11) shows very sim- 
ilar profiles; so that we decided to present the higher sen- 
sitivity observations of the isolated component. We also 
considered two different conditions in the model: (1) the 
N 2 H + abundance is constant throughout the core, so that 
the N 2 H + column density simply follows the dust, and (2) 
there is a "hole" (2000 AU in size) in the N 2 H+ distribu- 
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tion. The difference between the two synthetic profiles is 
not significant (only a 4% increase of the lincwidth to- 
ward the center, in the presence of the molecular hole, see 
Fig. |5J, after the convolution with a gaussian with a = 
y/ijfh + (T tu = 0.11 km s _1 , needed to match the intrinsic 
linewidth toward the center 4 , and thus only one profile is 
shown in Fig. El Also shown in the Fig. El are model 
profiles for the limiting case of no turbulent or thermal 
broadening (filled histograms). The thing to note is that 
the agreement with the data is mixed, in the sense that 
the predicted velocity gradient along the minor axis is ob- 
served but it is restricted to the south-west half of the 
axis. Along the major axis there is no clear gradient in 
the north-west half of the axis, as expected in absence of 
disk rotation, but a (0.06 km s _1 ) blue shift of the line is 
present toward south-east. 

The synthetic profiles become narrower as we move 
away from the central 40" of the disk-like cloud, given 
that the inward velocity in the adopted t% model reaches 
its maximum (0.12 km s _1 ) at radius 0.025 pc (or 37") be- 
fore rapidly dropping to at the cloud edge (see Fig. 2 of 
CB). We have analysed this line width variation to check 
its consistency with the observed trend shown in Fig. \§\ 
The solid curve in Figgis the CB prediction in the case of 
no central "hole" , whereas the dashed curve illustrates the 
effects on the line width of a molecular "hole" in the cen- 
tral 2000 AU (i.e. a region where all heavy elements have 
condensed onto dust grains). The line width of model lines 
has been calculated as the second moment of the velocity 
profile and it has been "normalized" to the value of the 
line width observed in the central position by convolving 
the purely-kinematic profiles with a Gaussian with a = 
0.11 km s _1 , which can be interpreted as the combina- 
tion of thermal broadening plus a constant turbulent field 
across the cloud. The comparison between the curves and 
the big dots (the binned data, see Sect 13. 5(1 suggests that 
our data are consistent with the CB model. In particular, 
the correspondence between the solid curve and the data 
indicates that the molecular hole is not present in L1521F, 
in agreement with the chemical analysis (see previous sec- 
tion) . On the other hand, the steeper Av-b relation found 
in L1544 (see Fig. 5 of Caselli et al. I2002a[> suggests that 
the molecular hole is likely to be present in that source, 
again in agreement with the chemical analysis (see e.g. 
Caselli et al. l2"U0l)l . 

As shown in Sect. 13.61 the kinematics of LI 52 IF is also 
characterized by complex motions which may be the re- 
sult of turbulence (see e.g. Burkert & Bodenheimer 120(F))) 
or accretion of material onto the core or a combination 

4 The linewidth toward the center and the observed Av-b 
trend can also be matched assuming a central infall speed 1.5 
times larger than the one calculated by CB in their model 
£3 and no extra broadening (e.g. turbulence) needs to be in- 
voked. However, this is not consistent with the CB model at 
time ts, and we do not further discuss this possibility. On the 
other hand, recent H2D + observations toward L1544 (Caselli 
et al. l2003il also imply larger central infall speed than predicted 
by CB-t 3 (van der Tak et al. EOPSll . 



of both. Local gradients have been determined in order 
to quantify these motions, and we found that the magni- 
tude of local gradient vectors tends to increase toward the 
core center. This suggests that turbulence (expected to 
dissipate more rapidly in denser environments), is proba- 
bly not the driving source of the observed velocity field. 
Moreover, unresolved substructure may further compli- 
cate the velocity field, as suggested by the N 2 H + clumps 
observed by Shinnaga et al. I|2004fl (see Fig. lllfl : with the 
exception of their clump "N4" , the kinematics of the other 
clumps (N1-N3) is consistent with the velocity field in- 
ferred from our local gradient maps (see also Sect. 13. lj) and 
the (marginal, see Fig. evidence for two velocity com- 
ponents in N2H + (3-2), resembling clumps Nl and N2). 

It is interesting to compare our data with predictions 
from the non-magnetic turbulent models of Ballesteros- 
Paredes et al. Ij2003|l . who derived velocity profiles for 
dense cloud cores. As an example, Fig. El shows the ve- 
locity cuts observed across L1521F. Within the half maxi- 
mum contour of N2H + (l-0), the largest velocity variation 
observed is 0.04 km s _1 on a scale of 0.027 pc, corre- 
sponding to a velocity gradient of ^ 1.5 km/s/pc. On the 
other hand, Ballesteros-Paredes et al. (2005|) found that 
the smallest velocity variation for their clump 13 at time 
tl (see their Fig. 9) is 0.3 km s _1 within 0.15 pc or > 2 
km/s/pc. Thus, current supersonic turbulent models pre- 
dict velocity gradients which are somewhat too large. We 
also note that the reversal in the velocity gradient direc- 
tion observed in L1521F (see bottom panel of Fig. I16f) is 
not present in the model examples shown in Ballesteros- 
Paredes et al. I|2003|) . but it is predicted by the turbulent 
models of Burkert & Bodenheimer ( 2000 ) . 

In the two starless cores LI 498 and L1517B, Tafalla 
et al. H2004II found a good correlation between the dis- 
tribution of CS and the distribution of the high velocity 
N2H + (l-0). The authors argue that the high velocity wing 
of the N2H + (l-0) lines comes from a gas shell that is be- 
ing accreted by the starless core; therefore it has not ex- 
perienced the high density in the core nucleus and hence 
its depletion factor is low. We repeated the same exper- 
iment in L1521F producing channel maps of the core in 
N2H + (l-0) but we did not find any strong deviation from 
the total intensity distribution. However, we did find that 
the N2H + (l-0) velocity pattern across the core is similar 
to the C 18 distribution. In fact, as shown in Fig. El 
N2H + (l-0) presents red-shifted velocities where CO is 
bright. The difference in magnitude of this effect here and 
in L1517B and L1498 could be due to a brighter core that 
saturates the emission from the high velocity wing. 

A final thing to note is that Fig. El together with 
Fig. ^2 (top figure), suggests the presence of a coherent ve- 
locity field resembling low-frequency spatial oscillations of 
the outer cloud layers around some equilibrium dynamical 
state as recently proposed by Lada et al. I|2003fl in the case 
of the starless core B68. However, L1521F is more mas- 
sive than B68 and is approaching the critical state for the 
onset of collapse, so that a situation of near-equilibrium 
for L1521F may be due to a balance between gravity and 
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a combination of magnetic and thermal forces. If this is 
the case, the normal modes for L1521F would be more 
complicated than in the purely thermal case of B68, since 
the spherical symmetry that is a fair approximation for 
a thermal-pressure supported cloud such as B68 will no 
longer be valid for a magnetically supported object (G. E. 
Ciolek, priv. comm.). In fact, we observe L1521F not to 
be spherically symmetric (see Tab.[3J|. 

4.3. Other possible interpretation 

Finally, in this section we discuss some more specula- 
tive interpretations of the evolutionary state of L1521F. 
Firstly, it seems that L1521F is likely to be in an earlier 
stage of evolution than L1544, which is probably colder 
(~ 7 K in the central ~ 2000 AU, see e.g. Zucconi et al. 
200 lj and more centrally concentrated (central densities 
~ 10 7 cm" 3 ; e.g. Evans et al. EMTjl than L1521F, and 
hence more depleted in the nucleus. For example, in the 
CB model, this implies that whereas the density profile of 
L1521F is consistent with the cloud model at time t 3 =2. 66 
Myr, L1544 is closer to t 4 =2.68 Myr, and thus the esti- 
mated age difference is roughly 20,000 yr. However, these 
age estimates should be taken with caution given that the 
two cores may have had a different formation and con- 
traction history and that, in general, cores, like stars, can 
differ in properties such as size, mass, and temperature 
regardless of their relative age. 

An alternative view can be based on the clumpy struc- 
ture observed in L1521F (Shinnaga et al. 2004), but not in 
L1544 (e.g. Williams et al. [TMfy . using BIMA. If the more 
complex kinematics in L1521F is due to the unresolved 
clumpy substructure, one may speculate that L1521F is 
close to the formation of a group of low mass protostars, 
whereas LI 544 is likely to form one or two stellar objects. 
Thus, the two cores may show chemical and physical prop- 
erties characteristic of the initial conditions of different 
modes of star formation in low mass cores. One should 
also note that L1521F resides in the middle of the main fil- 
ament of the Taurus Cloud, whereas L1544 is at the edges 
of the complex and hence the different environments of the 
two cores may be responsible for the different kinematics 
and chemical properties of apparently similar (in the dust 
continuum and C 18 emission) dense cores. 

5. Conclusions 

We have analysed physical and chemical properties of 
L1521F, a starless core in the Taurus Molecular Cloud, 
with characteristics similar to the pre-stellar core L1544. 
The main similarities and differences between the two 
cores are listed below: 

1. the dust emission distributions are similar, imply- 
ing a fairly closely matched density structure, with cen- 
tral densities of no ~ 10 6 cm -3 , the radius of the "flat" 
region ro = 20", and similar asymptotic power index a 
(see Sect. 13.21 and Tafalla et al. I2002[) . In particular, the 



aspect ratio is quite similar: 1.6 and 1.8 in L1521F and 
L1544, respectively. 

2. The line width decreases with distance from the 
cloud center (~ 0.3 km s _1 ) to 80" away (~ 0.25 km s , 
see Fig. in analogy with L1544, and consistent with 
the predictions of ambipolar diffusion models, although 
any gravity-driven contraction in 2-3D is expected to get 
localized line broadening. The particular model which best 
fits the data will be investigated in the future. 

3. The amount of CO freeze-out (integrated CO de- 
pletion factor /d = 15) is also comparable to L1544, as 
is the column density of N2H + toward the dust peak (~ 
1.5xl0 12 cm" 2 ). 

4. The deuterium fractionation toward the L1521F 
dust peak (i?dcut ~ 0.1) is however a factor ~ 2 smaller 
than in L1544, due to the (factor of 2) smaller column 
density of N2D 4 " toward L1521F. This can be understood 
if L1521F is less chemically evolved than L1544, with a 
smaller (r < 2000 AU) central molecular hole. 

5. Unlike in L1544, the velocity field in L1521F main- 
tains a complex structure even at the small scales traced 
by N 2 D+ and N 2 H+(3-2) (see Figs. HHEJ). This may be 
due to the presence of clumps in the central few thou- 
sand AU (as deduced by the interferometric observations 
of Shinnaga et al. 2004), but could also be caused by nor- 
mal mode pulsations, as in the case of B68 studied by 
Lada et al. (2003J). The ambipolar diffusion model with 
infall of Ciolek & Basu (2000) has problems in reproduc- 
ing the whole velocity field observed across the core. This 
may be due to the fact that part of the observed bulk mo- 
tions result from residual core contraction, as suggested by 
Tafalla et al. pMjl in their study of L1517B and L1498, 
thus preventing a clear determination of the velocity field 
within the core nucleus. 

6. The line profiles in L1521F show asymmetric struc- 
ture, although the two peaks are not well separated as in 
L1544. This is consistent with smaller (factor of 1.5) infall 
velocities in the central few thousand AU of L1521F. 

In any case, the large central density (~ 10 6 cm -3 ) 
and the evidence of central infall (broader line widths to- 
ward the center and central infall asymmetry in N2H + (1- 
0)) indicate that L1521F is another starless core on the 
verge of star formation, or a pre-stellar core. Although 
a study of a more complete sample is needed, assuming 
that L1544 and L1521F are the only two cores in Taurus 
close to dynamical collapse, and given that the total num- 
ber of starless cores in Taurus is 44 (Onishi et al. 2002), 
we can argue that this process of contraction towards the 
"critical" stage, or the "L1544-phase", may last about five 
percent of the core lifetime. 

A more detailed analysis of N2H + line profiles will be 
presented in a future paper, where the observed chemical 
abundances will be introduced in a non spherically sym- 
metric Monte Carlo radiative transfer code. This is needed 
to understand the nature of the N2H + (l-0) line asymme- 
try, which may be caused by self-absorption plus infall, or 
to infall plus a molecular hole, or to the relative motion of 
different clumps, or to a mixture of the above possibilities. 
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Fig. 4. Enlargement of the N 2 H+(3-2) spectrum of Fig. [31 
toward offset (-10,10), which shows tentative evidence for 
a double-peaked feature. The top panel show the line pro- 
file (grey histogram) and the Ms fit (black curve) assuming 
two velocity components; the bottom panel show the hfs 
fit assuming one velocity component. The velocities of the 
two components concide with the LSR velocities of clumps 
Nl and N2 observed by Shinnaga et al. ( 2004]) in their in- 
terferometric observations, see Fig. for reference. 

Fig. 5. Circularly averaged dust emission fitted with 
the expression no/(l+(r/ro) a ) convolved with a 
10.5"gaussian. The parameters that best fit the L1521F 
data are tiq = 10 6 cm" 3 , tq = 20"and a = 2. Dots with 
error bars are binned data and the dashed curve shows 
the observational beam. 



Fig. 1. Dust emission from L1521F. Levels are 30, 55 and 
80 mJy/beam. Reference position is RA: 04:28:39.8 DEC: 
26:51:35 (J2000). The dashed ellipse best fits the core 
structure with a 2D gaussian. The dotted ellipse is the 
result of a 2D-gaussian fit to the whole map, including 
the more extended emission. The black rectangle shows 
the area mapped in line observations (see Fig. 01. 



Fig. 6. CO depletion factor (full contours) against dust 
emission smoothed to a 22"beam (grey scale), /d con- 
tours range between 6 and 15, in steps of 3. The maximum 
value of /d is 18 and is located at offset (0,9). Note the 
good correspondence between the distribution of N 2 D + (2- 
1) (dashed contours) and the CO depletion map. 



Fig. 2. L1521F integrated intensity maps in the ob- 
served molecular transitions. Contour levels are 45, 70, 
95% of the relative peak in each map (whose val- 
ues are: 5.9, 0.74, 1.1, 0.33, 2.5 and 2.3 K km s _1 in 
N 2 H+(l-0), N 2 H+(3-2),N 2 D+(2-l),N 2 D+(3-2),C 18 O(l-0) 
and C 18 0(2-l) respectively). The N 2 H+(3-2) data were 
smoothed to 26"resolution to increase S/N and help 
the comparison with N 2 H + (l-0); for the same reasons 
N 2 D+(3-2) was smoothed to 16" (the HPBW of the 30m 
antenna at the frequency of the the N 2 D + (2-l) line). We 
overlaid on each map the ellipse which best fits the dust 
emission structure (in dashed line) and the 80 mJy/beam 
contour (dotted line). 



Fig. 3. Spectra of N 2 H+(l-0), N 2 H+(3-2), N 2 D+(2-l), 
N 2 D+(3-2), C 18 O(l-0), and C 18 0(2-l) towards the dust 
emission peak in L1521F. A fit of the hyperfine pattern 
(or in the case of C 18 a gauss fit) of the lines assuming 
identical excitation temperatures for all satellites is also 
shown (dashed) for comparison. 



Fig. 7. Integrated CO depletion factor /d against H 2 col- 
umn density as derived from the data in Fig.|SJ The dashed 
curve is the linear least fit to the data, whereas the solid 
curve is the /d vs. iV(H 2 ) relation found in the simple 
model described in Sect. I4.1T21 

Fig. 8. [N 2 D + ]/[N 2 H + ] in the central region as estimated 
from CTEX calculations (see Sect. \'6A\ . Levels are 0.066, 
0.078, and 0.091. The dashed contour is the small ellipse 
of Fig. El 



Fig. 9. N 2 H + (l-0) intrinsic line-width vs. impact param- 
eter. The big dots are averages of the data points within 
15" bins. The solid curve is the Aw — b relation derived 
from the velocity field predicted by the ts model of Ciolck 
& Basu (2000) (sec Caselli et al. 2002a), where a constant 
turbulent field across the cloud has been included, and the 
dashed curve is from the same model but accounting for 
the presence of a central molecular "hole" of 2000 AU in 
size (see text). 
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Fig. 10. Isolated component profile along the P.A. = 
45°axis passing through (-10,0) (left panel) and -45°axis 
passing through (-10,0) (right panel). This figure shows 
how the linewidth increases going from the edge of the 
core to the center. Fits are gaussian fits to this compo- 
nent. 

Fig. 11. Velocity gradient vectors in L1521F derived from 
the two N2H + maps. The integrated intensity of N2H + (1- 
0) (top) and N 2 H + (3-2) (bottom) is shown as grey scale. 
The thick contour marks the integrated intensity half- 
maximum contour. The arrows across the map are "local" 
velocity gradients obtained in a 3x3-point grid of posi- 
tions centered on the corresponding grid, indicating the 
magnitude and the direction. The arrow in the bottom 
right of the figure is the "total" velocity gradient derived 
from the whole set of available positions with A/ a a > 10 
(A being the integrated intensity). Note the large varia- 
tion in magnitude and direction of the local gradients. The 
white circles in the top figure locate the positions of the 
four N2H + clumps detected by Shinnaga et al. (2004). 

Fig. 12. Same as Fig.llllfor N2D + data. A similar pattern 
is present. 



Fig. 13. Deuterium fractionation (-Rdcut = 
A^(N 2 D+)/A r (N2D+)) as a function of integrated CO 
depletion factor. Filled squares are data from this paper 
and the filled circle with error bar is our result for the 
dust peak of L1544 (Caselli et al. l2002b» . Thin full curves 
are predictions from a simple chemical model using dif- 
ferent rate coefficients for the proton-deuteron exchange 
reaction Q: the "standard" rate of 1.5xl0~ 9 cm 3 s _1 
and the " GHR" rate of 3.5X10" 10 cm 3 s -1 , recently 
determined by Gerlich et al. ( 2002). Dashed curves refer 
to the same models but with larger value of a m ; n , the 
minimum radius of dust grains in the adopted MRN 
grain-size distribution (see Sect. 14. lTTfl . The dotted curve 
is the best fit to the data, found if fcjjD = 7.5xl0~ 10 cm 3 
s . The thick curve is the result of a more comprehen- 
sive chemical model which takes into account the density 
structure of the core (see Sect. I4.l"2)l . Typical l-cr error 
bars are shown. To reproduce the i?deut value observed 
in L1544, N 2 freeze-out and multiply deuterated forms 
of H^" have to be included in the model (dash-dotted 
curve) . 



Fig. 14. Fractional abundances, n(i) / n(H2) , as a function 
of radius in L1521F. The symbol M + refers to metals (see 
text). The adopted density profile is the one described in 
Sect . 13.21 This model reproduces the observed column den- 
sities of CO, N 2 H+, and N 2 D+. The N 2 H+ abundance de- 
creases towards the center by a factor of about 2, whereas 
N 2 D + increases by a factor of ^7 from core edge to core 
center. Note that the H 2 D + abundance is predicted to be 
~3xl0 -10 at the core center, a factor of 3 smaller than 
observed in L1544. In analogy with L1544, the HCO + 
and DCO + abundance profiles present a steep drop within 
about 5000 AU. The right axis refers to the CO depletion 
factor (see the curve labelled "Fb")- The set of param- 
eters used to obtain this model are: E'd(CO) = 1210 K, 
£ D (N 2 ) = 800 K, £b(0) = 750 K, C = 1.3xl0- 17 a -1 , 
Omin = 5xl0~ 6 cm, and a;(N 2 ) = 2xl0~ 5 , where En(i) 
is the binding energy of species i and £ is the cosmic-ray 
ionization rate (see Caselli et al. 2002b for details of the 
model). 

Fig. 15. Observed line profiles of the F 1 F = 1^1 2 
hyperfine component of N2H + (l-0) (empty histograms) 
along the minor (left panels) and major (right panels) axis 
of L1521F compared with predictions based on the Ciolek 
& Basu 2000; velocity field at time t = (full curves), 
convolved with a a = 0.11 km s" 1 gaussian to reproduce 
the observed intrinsic linewidth of the central spectrum 
(see Fig. [2J . The filled histograms are model line profiles 
computed assuming no turbulent or thermal broadening. 
The core center is defined by the dust peak, located at 
offset (-8,+2). 

Fig. 16. Intensity (thin lines) and velocity (thick lines) 
of N2H + (l-0) along four different cuts; in the top panel 
continuous lines refer to the minor axis (PA 115°) while 
the dashed line refers to a cut at PA 70°; in the bottom 
panel, the continuous lines show the behaviour along the 
major axis (PA 25°) and the dashed lines represent the 
other bisector (PA 160°). 

Fig. 17. C 18 O(l-0) integrated intensity map (greyscale; 
wedge units are K km s _1 ) overlaid by the N 2 H+(l-0) 
velocity map (contours). N2H+ velocity levels are from 
6.37km s^ 1 to 6.55km s~ 1 spaced by 0.02km s _1 ; higher 
velocity contours are represented as solid lines. 
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